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Abstract 

The diffusive characteristics of two upwind schemes, multi-dimensional fluctuation split- 
ting and locally one-dimensional finite volume, are compared for scalar advection-diffusion 
problems. Algorithms for the two schemes are developed for node-based data represen- 
tation on median-dual meshes associated with unstructured triangulations in two spatial 
dimensions. Four model equations are considered: linear advection, non-linear advection, 
diffusion, and advection-diffusion. Modular coding is employed to isolate the effects of the 
two approaches for upwind flux evaluation, allowing for head-to-head accuracy and efficiency 
comparisons. Both the stability of compressive limiters and the amount of artificial diffusion 
generated by the schemes is found to be grid-orientation dependent, with the fluctuation 
splitting scheme producing less artificial diffusion than the finite volume scheme. Conver- 
gence rates are compared for the combined advection-diffusion problem, with a speedup of 
2.5 seen for fluctuation splitting versus finite volume when solved on the same mesh. How- 
ever, accurate solutions to problems with small diffusion coefficients can be achieved on 
coarser meshes using fluctuation splitting rather than finite volume, so that when comparing 
convergence rates to reach a given accuracy, fluctuation splitting shows a speedup of 29 over 
finite volume. 


Nomenclature 


A, B Flux Jacobians 

c Nodal update coefficients 

F Convective flux vector, F = F[x, y, u) 
h Mesh spacing 

i,j Cartesian unit vectors 

J -1 Inverse Jacobian of the coordinate transformation 
i Edge length 

M Averaging function for limiters 

n Outward unit normal to control cell 

n Length-scaled inward normal 

Q Fluctuation ratio 

f Distance vector from node to face 

Si Median-dual area about node i 

t Time 

u Dependent variable 

x, y Cartesian coordinates 

a, 0 Curvilinear advection speeds 

T Boundary of control cell 

7 Limiter bound 

e Small parameter for van Albada limiter 

d Finite element linear shape function 

A Advection velocity vector 
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V 

Diffusion coefficient 


Curvilinear coordinates 

T 

Timestep 

ft, ft 

Artificial dissipation — fluctuation splitting 

$ 

Artificial dissipation — finite volume 

<t> 

Element advective fluctuation 

ft,4> n 

Fluctuation components 

ft' ft' 

Limited fluctuations 

4>v 

Viscous fluctuation 

V 

Limiter function 

n 

Area of control cell 

V 

Gradient operator 


Introduction 

Upwind discretizations for advection equations typically introduce artificial numerical 
dissipation into the solution. When combined advection-diffusion problems are considered, 
this dissipation introduced in the discretization of the advection terms should be less than 
the true physical diffusion. To this end the diffusive characteristics of upwind schemes are 
investigated and their performance in resolving solutions to advection-diffusion problems 
with small diffusion coefficients is analyzed. 

Two node-based, median-dual methods for modeling convective fluxes are considered. 
The first, is a traditional locally one-dimensional approximate Riemann solver finite volume 
(F\ ) scheme . 1 Locally one-dimensional schemes applied on multidimensional domains are 
known to introduce excess dissipation when discontinuities are not aligned with the mesh . 2 

The second method is the NNL 3 fluctuation splitting (FS) scheme, also referred to in 
the literature as a residual distribution scheme. FS has a more-compact stencil than FV for 
second-order formulations and exhibits ‘‘zero cross-diffusion” t in a grid-aligned condition. 
Both of these attributes should lead to less introduced dissipation as compared with FV. 

The sensitivity of FS to grid orientation and resulting production of cross-diffusion is 
investigated in the present report. The use of compressive limiter functions is also tested 
with both algorithms. Local timesteps based on positivity arguments are implemented for 
both first- and second-order discretizations of the implicit matrix. 

Formulation of FS schemes for diffusion problems is a recent research area . 4,5 The present 
study seeks to quantify the relative merits of using i low-diffusion advection operator to 
resolve advection-diffusion problems with small diffusion coefficients. Lessons learned on 
these problems will guide the development of computer ^odes for solving compressible viscous 
fluid dynamic problems. A similar approach for central difference schemes with explicit 
numerical dissipation has recently been taken by Efraimsson . 6 

4 “Zero cross diffusion” refers to the practice of adding artificial diffusion terms in the streamwise direction 
only, as opposed to adding artificial dissipation in both the streamwise and cross-stream directions. 
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Governing Equations 

The non-linear advection-diffusion equation, 

u t + VF = V • (i/Vm) (1) 

is cast as a hyperbolic conservation law, to which steady-state solutions are sought. 

Finite Volume 

In FV form, using the divergence theorem Eqn. 1 becomes, 

I u t dQ = - j {F - i/Vii) • fi dT (2) 

where Q is the median dual about node i and T is the boundary of Q. Using mass lumping 
to the nodes, similar to an explicit finite element treatment, 8 the temporal evolution is 
evaluated on a time-invariant mesh as, 

j u t dQ = 5 , —I — > — («| +T — u'i) (3) 


The discretization of the convective flux, F, is performed using Barth’s implementation 1 
of the upwind, locally one-dimensional, approximate Riemann solver of Roe. 9 


J) F ■ h dV — > 'y ] — ^ F m + Fgut'j 

' faces 


n 


$ 


AT 


( 4 ) 


where the artificial dissipation provides the upwinding, 

$ = — T I$iiy\(7i ou i Uj n ) (5) 

with h = h x i + hyj. Out and in refer to states on the outside and inside of Q at the face. A 
and B are the flux Jacobians, 


dF(l) 2) 

^ “ du ’ ” du 

and (A, B) represent their conservative linearizations at the cell face. 9 

Piecewise linear reconstruction from the nodal unknowns to the cell faces as, 


( 6 ) 


Uface = Uj + i’Vu ■ T (7) 

provides second-order spatial accuracy in smoothly-varying regions of the solution. Median- 
dual gradients of the dependent variable, Vu, are obtained from the unweighted least squares 
procedure outlined by Barth. Following Bruner and Walters, 14 the limiter is supplied an 
argument equal to half the argument Barth uses, namely, 

( yminjmax _ 

— — — . , 1 

2(\7 U * fA rn ^/max 
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where u jmn / max i s the minimum (resp. maximum) of and all distance-one neighbors. The 
most restrictive limiting from choosing the minimum or maximum is used. 

In casting the limiter argument in this form, Bruner equates the Barth limiter with 
Superbee, for a limiter argument less than or equal to one. The present authors incorrectly 
identified the Barth limiting with the non-symmetric Chakravarthv and Osher 13 limiter in 
Ref. 7. The Barth limiting is non-symmetric, but takes the form, 

( 0 | < 0 

1 2; if «<?<! (9) 

l 1 

for the limiter cast as Eqn. 8. 

Tw'o methods for evaluating the diffusion term are incorporated into FV. The more 
compact of the two, the finite element discretization, is discussed in the following section. 
The less-compact, diffusion formula is obtained by discretizing the last term of Eqn. 2, in a 
manner similar to Eqn. 4, 

f uVu • ndT -> ^ (V« jn 4- Vu^) ■ h A Y (10) 

Jv faces 2 

The diffusion coefficient is averaged over the length of the face. The gradients from Eqn. 7 
are not limited before averaging at the control-volume faces in Eqn. 10, as suggested by 
Anderson and Bonhaus. 10 



Fluctuation Splitting 

The NNL FS scheme is presented as a slight re-interpretation of the work of Sidilkover and 
Roe. 3 The current interpretation is as a volume integral over triangular elements, without 
recourse to the divergence theorem. The discretized equations, however, are identical. 
Integrating Eqn. 1 over an element, where Q is now the area of the triangular element, 

u t dQ = - I V • F dD + I V • (1 sVu) dil (11) 

./ n Jn Ju 

For linear variation of the dependent variable over the element, the temporal evolution is, 



Qu t = 


d . 

■g(«l« T U ‘2, +tl 3( ) 


( 12 ) 


where u 2 , and 113 correspond to the three nodes defining element 0. 

Defining local curvilinear coordinates, £ and r/, parallel to sides 12 and 23, respectively 
(Fig. 1), the divergence of the convective flux can be written, 


v • F = fi" + Ff = Af (n 2 . F f - ft, . F„) (13) 

Defining the scaled inward normal, n = —hn, where h is a mesh edge length, the divergence 
(Eqn. 13) becomes, 


V • F = 


2D 


(-/ii 2 n 2 • + /t 23 n! • 


(14) 
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If F is linear or quadratic in u, then for a linear variation of u over the element. 


V • F dil = aA‘2iu + pA 32 u 


J n 

where the difference operator is defined A- 2 \U = u 2 — U\ and the advection speeds are, 

° + P = ^(n lx i + n li( B) 

A and B are now the conservative linearizations over the triangular element . 11 
The advective fluctuation can be defined, 


(15) 

(16) 


The fluctuation can be split, 


L 


0 = - / V FdQ 


(17) 


0 = 4 F + 4> v 


(18) 


wliere, 

0* = -aA 2 iu, 4> ] = -PA 32 u (19) 

Following Sidilkover 12 the fluctuation is limited to achieve a second-order scheme, 

pp = + P v i>{Q) = (ft ^1 - ( 20 ) 

<t> v ' = P" - 0V(Q) = P" (1 - P(Q)) (21) 

with, 

Q = (22) 

Upwinding is achieved through the introduction of the artificial dissipation terms, 

0 ^ = sign(a)0^ , p n = sign (p)p T> (23) 


Combining Eqn. 12 with a distribution scheme for Eqn. 17 and summing over all elements, 
the contributions to nodal time derivatives can be written in tin 1 form, 

Siu u <- \(P V -&) + COE 

S 2 u 2t <- \(P*' + P*) + \(P V ' -pv + coe 

S 3 u 3 , <- l -(r' +r) + COE (24) 
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or in a more compact form, 


Sun , <- ^ [*(3 - i)(<t> v + (-1 Yft) + (-4 + 5 i- z 2 )(<^‘ - (-1 )>")] + COE (25) 

where COE stands for contributions from other elements containing these nodes. 

A finite element treatment, similar to Tomaich, 4 is employed to obtain the diffusive 
fluctuation, 

(j>v= / V(*/Vu)dH ( 26 ) 

Jn 

Assuming piecewise-linear data and an element-averaged diffusion coefficient leads to a dif- 
fusive fluctuation of zero for the triangular element. Introducing the linear nodal shape 
functions i?j, such that — E the elemental diffusive fluctuation can be expressed 

(pv = Yli = 1 <t>vi = 0* where 

<t> Vi — I t?fV*(PVM)rffi ( 27 ) 

Jn 

Integrating by parts, 


<t>vi = j> tfjt'V u ■ h dF — j i/Vn • Vtf* dVl (28) 

The boundary integral in Eqn. 28 will cancel on summing contributions for interior nodes. 
The remaining volume integral can be evaluated analytically, 


3 

K = ~2 Vu ' n «+i = -Tq “t n j+i • n i+i (29) 

j~ 1 

Distributing this diffusive fluctuation to the nodes and keeping only the larger of the physical 
or artificial dissipation leads to the update formula, 


~ <t> V 

oit/q < — — + max 


+ <j> v 

>> 2 ^ 2 , < 1- max 


c (r 

*^3«3 , <- y + max I y 


0 1 \ 

y,<M +CO£ 

dx ( 2 ’ J 

,<0 +COE 


+ COE 


(30) 


Boundary Conditions 

Explicit Dirichlet inflow boundary conditions are employed. Advective outflow bound- 
aries are treated for free convection through the boundary nodes, allowing boundary nodes 
to be handled in the same manner as interior nodes. For the diffusion terms a Neumann 
outflow boundary is applied with zero gradient, achie/ed by setting the boundary integral 
in Eqn. 28 to zero. 
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Limiter functions 

Minmod, van Albada, 15 Superbee, and 7 16 symmetric limiters are utilized for FV (Eqn. 7) 
and FS (Eqns. 20 and 21) in the form of symmetric averaging functions related to the limiter 
as, 

Q ip ^ = M (p, q) = M (q, p ) = p J 


The van Albada averaging function is, 


(pq + g 2 )(p + g) 

p 2 + q 2 + 2e 2 


where the small parameter £ varies like £ 2 ~ Ax 3 , and serves to reduce the limiting in smooth 
regions. 

The averaging function for the 7 limiter, of which the Minmod (7 = 1 ) and Superbee 
(7 = 2) are special cases, is, 


M(p,q) 


0 pq < 0 

7 P 7|p|<kl 
< 9 if |p| <kl < 7|pI 
P \q\<\p\< l\q\ 

, 79 7kl<|p| 


( 31 ) 


Timestep 

Both schemes are formulated either as Gauss-Seidel time-relaxation or forward Euler 
time-evolution algorithms. 

The nodal updates for the discrete system can be formed as a sum of contributions from 
all nodes. 


u\ +T = ^2 c j u j = CiUi + 22 C j u j ( 32 ) 

j 

For positivity 17 each of the coefficients in Eqn. 32 must be non-negative. 

Advective Timestep restriction 

In the FV context the nodal update (Eqn. 32) can be rearranged into the form of Eqn. 3, 

^0h- +T - u l) = y (Ci - l)«j + y 22 c 3 u i ( 33 ) 

For the upwind, edge-based algorithm considered here, each 2 c j will be related to a positive- 
definite coefficient equal to zero for outflowing faces and related to the wavespeed for in- 
flowing faces, yielding the restriction r > 0 on the timestep. The remaining term can be 
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expressed 


f(«-l) = -5> (34) 

k about i 

where the c k coefficients are also positive-definite, either zero for inflowing faces or related 
to the wavespeed for outflowing faces. Rearranging and imposing the positivity constraint, 
<i > 0, yields the timestep restriction, 

1 - j X/* = c « ^ 0 (35) 

k about i 


T - v T ( 3 <i) 

Xsk about i K 

For FS, the nodal updates are assembled from Eqn. 24 as, 

7 ( U l +T ~ U D = X c t(*b - u ‘) (37) 

In this case the cj coefficients are formed as contributions from the fluctuations in the 
triangles to both the left and the right of mesh edge rj. The positivity restriction on r is 
found to have a similar form as for finite volume (Eqn. 36), 


Sj 

c j 


(38) 


Local time-stepping based on positivity is showm to yield stable, yet non-converging, so- 
lutions in some second-order cases (see Results section). Robust convergence is obtained by 
using the first-order c s in Eqns. 36 and 38, even for second-order-accurate spatial discretiza- 
tions. This is equivalent to the common practice of using a first-order Jacobian discretization 
in an time-implicit scheme. 


Diffusive Timestep Restriction 

Unfortunately, the finite element formulation for the diffusive terms (Eqn. 29) cannot be 
guaranteed to preserve local positivity on obtuse triangles (see Barth 1 ). Considering only 
the contributions from the current node, the coefficient for the diffusion term can be w'ritten, 


u 


t + T 




(39) 


The appropriate edge length is the side of the element opposite the current node. The 
resulting timestep restriction is, 


T < 


s t 


i/i 2 


St 4$i 

In a similar manner the timestep restriction from Eqn 10 is, 

5, 


T < 


Ei 




(40) 


(41) 
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Results 


Linear Advection 

The linear advection equation is obtained from Eqn. 1 by setting v = 0 and F = Xu, 
yielding, 

u t + V • (Au) = 0 (42) 

A divergence-less advection velocity is considered, such that V • A = 0. Equation 42 can then 
be written, 

Ut + A • Vu = 0 (43) 

Uniform Advection 

Uniform advection of the Heavyside function at —45 degrees, A — (1,-1), on a cut- 
cartesian mesh is shown for first-order FS, second-order FS, and second-order FV in Figs. 2 -4, 
respectively. The mesh is shown as the dashed background, and equally-spaced contours vary 
on [0,1], the minimum and maximum solution values. The spread of the contour lines with 
spatial evolution is indicative of the amount of dissipation introduced into the solution by 
the discretization of the convective terms. 

Second-order FS is seen to be greatly superior to first-order, as expected, reproducing the 
exact solution in this case with no introduced dissipation. Also, FS represents a significant 
reduction in numerical diffusion versus the corresponding FV scheme, with both results 
employing the Minmod limiter. 

However, the “zero cross-diffusion” results of Fig. 3 with FS are misleading. In Fig. 5 
the advection velocity has been rotated counter clockwise by 90 degrees on the same grid. 
Clearly, the artificial dissipation introduced by the FS scheme has been increased. 

The corresponding FV solution is shown in Fig. 6. While the change in contour spreading 
for the FV scheme between Figs. 4 and 6 is less dramatic than the change in spreading for 
the FS scheme in Figs. 3 and 5, the FS results still exhibit less diffusion than the FV results, 
comparing Figs. 5 and 6. 

Employing the compressive Superbee limiter with the FS scheme yields the results of 
Fig. 7. In this case the discontinuity is confined to a 2-3 cell stencil, and does not grow in 
space. Applying the Superbee limiter to FV cannot eliminate all artificial dissipation for 
this case, as is possible with FS. The FV results (not showui) corresponding to Fig. 7 spread 
the discontinuity over four cells by the outflow boundary, with a continually broadening 
trend. 

However, while it is possible to use the Superbee limiter with FS for this case, compressive 
limiters can be unstable on different grid orientations. For example, no degree of compression 
is stable for the case of Fig. 3. This potential for instability is related to global positivity, 
as discussed by Sidilkover and Roe. 3 

The effect of using a general unstructured grid is investigated in Figs. 8 and 9. The 
unstructured grid in this case was generated using Vgrid. 18,19 The FS solution exhibits 
less dissipation, but is not as smooth as the FV solution. While the FS scheme preserves 
contact discontinuities over larger spatial ranges than the FV scheme, FS does not appear 
to degenerate gracefully with regard to extreme coarsening of the unstructured mesh for 
this test case. This behavior could have negative implications for applications employing 
multigrid convergence acceleration. 
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Circular Advection 

Circular advection is achieved by setting A = (y, - r). A decaying sine- wave input profile 
is used. 


u(x , 0) = (e x sin nx) 2 

Results for the two schemes, using the Minmod limiter, are presented on the worse-case 
cut-cartesian mesh in Figs. 10 and 11. Again, the FS results are considerably less diffusive 
than the FV solution. 

The circular-advection problem is also applied on an unstructured mesh. The input 
profile for this case consists of both a top-hat function and a decaying sine wave, allowing 
comparisons between the schemes for both sharp discontinuities and smooth gradients. The 
input profile is, 

0.5 < x < 0 
0.6 < x < —0.5 
0.8 < x < -0.6 
-\<x< -0.8 

Results for this case are displayed in Fig. 12 for FS and Fig. 13 for FV, both using the 
Minmod limiter. FS performs significantly better at preserving the top-hat, distribution. 
FS also does a better job of maintaining the minimum and maximum values of the sine 
distribution, though both schemes do well on the smooth gradient portion of the sine wave. 

Non-linear Advection 

The non-linear advection equation is obtained from Eqn. 1 by setting F = (y, u) with 
v = 0. In non-conservative form the equation is written, 

u t + uu x + u y = 0 

A coalescing shock problem is considered, with an anti-symmetric input profile, 

«(-l ,y) = “(0 ,y) = 0 


(e 2x sin(27ra:)) — 


u{x , 0 ) = < 


0 

0.4 

0 


u(x, 0) — —2x — 1 on x = (—1,0) 

The exact solution to this problem contains symmetric expansion fans on the sides and a 
compression fan at the inflow that coalesces into a vertical shock at (x,y) - (-1, 1). 

The first mesh is cut-cartesian containing 26 x 2(1 nodes. The FS and FV solutions, 
both using the Minmod limiter, are presented in Figs. 14 and 15. Both algorithms exhibit 
the same grid dependence on the amount of artificial dissipation as seen before, with the 
left-half solutions having more diffusion than the righ , halves, due to the grid orientation. 
Both methods perform the same in the compression-f in region, coalescing into a shock to 
within the accuracy of the input-profile discretization. 

The shock is more sharply defined by FS than bv FV. Figure 14 has the correct shock 
speed, with nearly the entire gradient captured in one cell thickness. In contrast, Fig. 15 
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shows a slightly incorrect shock speed when using FV, as the shock progresses to the left 
beyond the coalescence point, even though the discretization is conservative. The incorrect 
shock speed results from a noil-symmetric distribution of the dependent variable to the left 
and right of the shock, caused by the excessive artificial diffusion generated on the grid- 
misaligned (left-hand) side. 

Contours of the absolute value of the error are presented in Figs. 16 and 17. Errors from 
both computed solutions show a lack of symmetry, again reflecting the grid dependence of 
the artificial diffusion terms. The error levels from FS are less than from FV. The shock 
curvature in the FV solution at the coalescing point is clearly visible in Fig. 17, resulting in 
significant downstream errors in the shock location as compared with the FS errors. 

This problem is repeated on a 25 x 25 mesh with symmetric diagonal cuts, favorably 
aligned with the advection directions. The FS and FV solutions, Figs. 18 and 19, are in 
good agreement. Plots of the absolute error contours, Figs. 20 and 21, show FS to be a 
little more accurate than FV for this case. 

The final mesh for this case is a truly unstructured triangulation containing 847 nodes 
and 1617 cells. The nodes are clustered to the outflow boundary, with a bias towards the 
left-hand side. The FS solution is presented in Fig. 22, showing very accurate and crisp 
shock resolution and good symmetry in the solution contours despite the mesh-clustering 
bias. In contrast, the FV solution in Fig. 23 has a more-diffuse shock and again an incorrect 
shock speed, with the outflow' shock offset to the left of x = — The FV solution is also 
somewhat less symmetric than the FS solution. 

Linear Diffusion 

Choosing F = 0, the heat-conduction equation is obtained from Eqn. 1, 

u t = V • (izVu) 

The test problem, a steady-state boundary value problem on a unit square, is taken from 
Tomaich. 4 The Dirichlet boundary values are, 

u{-l, y) = 0, u( 0, y ) = sin(7ry) 

u(x, 0) = 0, u(x, 1 ) = — sin(7r.r) 

The analytical solution oni = [—1, 0], y — [0, 1] is, 

u(x, y) = — ^ — fsinh(7r(x+l)) sin(7ry) + sinh(7ry) sin(7r (rr-h 1))] 

SUlll 7T 

Both diffusion discretizations, Eqns. 10 and 29, are compared on a 438-node unstructured 
mesh. Figures 24 and 25 plot the absolute value of the error in the converged solutions 
using Eqns. 10 and 29, respectively. A carpet plot of the solution, using the finite element 
formulation, is presented in Fig. 26. 

The finite element treatment is clearly more accurate, and is used to discretize the dif- 
fusion terms for both FV and FS in the following section. The average-gradient, results in 
Fig. 24 appear to exhibit a decoupling mode, similar to odd/even decoupling for structured 
meshes. 
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Linear Advection-Diffusion 

The final test^ case is a linear advection-diffusion problem of Smith and Hutton. 20 The 
flux function is F — Xu, with, 

A = (2y(l - x 2 ), -2x{] - y 2 )) 

The streamlines for this problem, while not truly circular, are similar in orientation to the 
circular advection problem. The inflow profile is, 

u(x, 0) = 1 + tanh(20x + 10) 

The diffusion coefficient is chosen to be a constant, v = 10 -3 . The domain is the unit square 
in the second quadrant. No closed-form solution is known to the authors for this problem. 

A sequence of five unstructured meshes is considered. The meshes have no preferred 
clustering or stretching and have nominal node-spacings of 0.1, 0.05, 0.025, 0.0125, and 
0.00625, labeled as Meshes A- E, respectively. The number of nodes for each mesh, along 
with the solution times for both FS and FV on a 195 MHz SGI R10000 CPU are listed in 
Table 1. 

Z/ 2 -norrns of the artificial and physical viscosities computed using both FS and FV are 
presented for each mesh in Table 2. Notice that the norm of the artificial dissipation 
for both FV and FS drops lower than the norm of the physical dissipation on Meshes D 
and E. Since the algorithms select only the larger of the physical or artificial dissipation 
(Eqn. 30), Table 2 suggests both schemes are grid resolved on Mesh D. However, the norm 
of the physical dissipation is smaller for FV than FS on each mesh A D. The physical 
viscosity is driven by the solution curvature, suggesting FS maintains the solution profile 
sharper than FV on the coarser meshes. A comparison of outflow profiles will soon verify 
this interpretation. 

Further evidence of a grid-resolved FS solution is seen in Figs. 27 and 28. The FS 
solution on Mesh E at the outflow boundary is presei ted along with the inflow profile and 
the corresponding pure-advection (u = 0) FS solution in Fig. 27. The pure-advection 
solution is seen to replicate the inflow profile, with a clear separation from the diffused, 
v = 10 \ solution. Plotting only the FS results with respect to grid refinement, Fig. 28 
shows a convergence of the outflow profile by Mesh C for FS. 

The accuracy of FS and FV are compared in Fig. 29, where the outflow solutions from FS 
and FV are plotted for Meshes C and E. Taking the grid-converged FS Mesh-E solution to 
be the '"truth” solution, it is clear that FS reaches tin grid converged solution on a coarser 
mesh than FV. 

Computational efficiencies of the two algorithms aie compared in Fig. 30, where the L 2 - 
norm of the residual is plotted versus CPU time for the fine-mesh FS and FV solutions, along 
with the FS convergence history on Mesh D. The Mesh-E FS solution converges in 760 sec. 
The corresponding FV solution takes 2.5 times longer than FS, due, in part, to the need 
to reconstruct gradient information at each node with FV for second-order spatial accuracy. 
However, considering the solution time to reach a given accuracy, it is more reasonable to 
compare the FS solution time on Mesh D to the fines -mesh FV solution. The FS Mesh-D 
solution took only 64 sec, a factor of 29 times less than FV on Mesh E, and still shows better 
accuracy than the fine-mesh FV solution. 
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An even greater speedup is seen with FS in conjunction with the van Albada limiter, 
where now the Mesh-B solution over-plots the curve from the finest grid, shown in Fig. 31. 
The corresponding FV result using the van Albada limiter on Mesh B is included, and 
clearly falls short of the FS accuracy. The FV case was repeated with the highly-compressive 
Superbee limiter with little improvement in accuracy. The solution time for FS on Mesh B 
is about one second, yielding a speedup factor of 2 3 orders of magnitude over FV. 

The final set of results addresses convergence issues while pushing the positivity limits. 
Figure 32 compares two convergence histories for the second-order FS on Mesh B. The 
non-converging, though stable, convergence history is the result of using strict positivity ar- 
guments to set the timestep (Eqn. 38). The resulting solution is bounded and approximately 
correct but oscillatory. Limiter “ringing” is considered to be a contributor to this behavior, 
and the higher-order discretization for the implicit matrix could be reducing the diagonal 
dominance, and hence stability, of the Gauss-Siedel iteration. 

Full convergence is achieved by using first-order positivity coefficients, which are not 
dependent on the limiters. The resulting local timesteps will not be as large as true second- 
order positivity would allow, but appear to be more robust. 

Summary of Results 

Fluctuation splitting and finite volume schemes are compared in detail as applied to scalar 
advection, diffusion, and advection-diffusion problems. The fluctuation splitting scheme is 
seen to introduce less artificial dissipation while treating advection terms, allowing for more 
accurate resolution of weakly dissipative advection-diffusion problems. The ability to resolve 
solutions to these problems on coarser meshes makes the fluctuation splitting scheme the 
preferred choice over finite volume. 

Linear advection test problems are utilized to investigate the dependence of artificial 
diffusion production on grid orientation. Both fluctuation splitting and finite volume are 
shown to exhibit grid dependencies, but with fluctuation splitting producing less artificial 
dissipation on all grids considered. 

A non-linear coalescing shock problem further explores grid dependencies as cases are 
constructed that result in incorrect shock speeds for finite volume. Fluctuation splitting 
shows correct shock speeds for all grids and provides tighter shock capturing than finite 
volume. 

An advection-diffusion problem with small physical dissipation (diffusion coefficient of 
10 -3 ) is considered where the reduction in artificial dissipation with fluctuation splitting 
results in a significant accuracy improvement over finite volume. Convergence times are 
compared, showing a speedup of 2.5 for fluctuation splitting over finite volume on identical 
grids, using a point Gauss-Seidel relaxation. However, a grid convergence study shows 
fluctuation splitting has better resolution of the solution on a coarser mesh than finite volume 
does on finer meshes, resulting is a speedup of 29 for fluctuation splitting over finite volume. 

Based upon these significantly reduced solution times for solving model problems, as com- 
pared to the current state-of-the-art finite volume method, fluctuation splitting is considered 
a worthwhile scheme to pursue for modeling fluid dynamic problems. 
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Table 1 Grids and solution times for ad vect ion-d iffusion problem. 

CPU seconds 
Mesh Nodes FS FV 
A 134 < 1 <1 

B 495 1 1 

C 1,928 5 8 

D 7,529 64 145 

E 28,915 760 1880 
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Table 2 L 2 -norms (xlO°) of artificial and physical viscosities for ad vect ion-diffusion 
problem. 


mu 

(art.) 

FS 

IM 

(phvs.) 

Mesh 

(art.) 

FV 

\\4>v\\2 

(phys.) 

1274 

215 

A 

1918 

190 

597 

265 

B 

640 

176 

192 

161 

C 

144 

119 

54 

76 

D 

46 

66 

13 

36 

E 

18 

36 
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Figure 2 First-order fluctuation splitting, uniform advection. 
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Figure 3 Second-order fluctuation splitting, uniform advection. 
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Figure 4 Second-order finite volume, uniform advection 
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Figure 6 Second-order finite volume, uniform advection 
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Figure 7 Second-order fluctuation splitting with compressive limiter. 
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Figure 8 Fluctuation splitting on unstructured mesh 




Figure 9 Finite volume on unstructured mesh. 




Figure 10 Fluctuation splitting, circular advection. 
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Figure 13 Finite volume on unstructured mesh, circular advection 
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Figure 15 Finite volume, Burgers equation. 






Figure 18 Fluctuation splitting, Burgers equation, symmetric mesh 
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Figure 22 Fluctuation splitting, Burgers equation, unstructured mesh 
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Figure 23 Finite volume, Burgers equation, unstructured mesh 
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Figure 24 Pure-diffusion problem 



X 


, diffusion terms from Eqn. 10. 




Figure 25 Pure-diffusion problem error, diffusion terms from Eqn. 29 





Figure 26 Heat equation solution using finite element formulation. Contour increment 
is 0.1. 
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Figure 27 Fluctuation splitting profiles on finest mesh, advection-diffusion problem. 
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Figure 29 Fluctuation splitting and finite volume for advection-diffusion problem. 




Figure 30 Convergence histories for advection-diffusion problem. 
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Figure 31 Advection-diffusion results using van Albada limiter. 
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